/*******************************************************************************
 2HDMC - two-Higgs-doublet model calculator
 THDMCScan: Based on Demo program

 http://2hdmc.hepforge.org
*******************************************************************************/
#include "Constraints.h"
#include "DecayTable.h"
#include "HBHS.h"
#include "SM.h"
#include "THDM.h"
#include <iostream>
#include <iomanip>
#include <string>
#include <sstream>
#include <fstream>
#include <cmath>

using namespace std;

int main(int argc, char *argv[]) {

  if ( argc != 15)
  {
	 printf("THDMCScan usage:\n");
	 printf("./THDMCScan <tag> <BitWriteLHA> <mh> <mH>  <mA> <tanb> <mC> <sinba> <l6> <l7> <m12_2> <yukawa-type>");
  }

  double mh_in    = (double)atof(argv[1]);
  double mH_in    = (double)atof(argv[2]);
  double mA_in    = (double)atof(argv[3]);
  double mC_in    = (double)atof(argv[4]);
  double sba_in   = (double)atof(argv[5]);
  double tanb_in  = (double)atof(argv[6]);
  double l6_in    = (double)atof(argv[7]);
  double l7_in    = (double)atof(argv[8]);
  double m12_2_in = (double)atof(argv[9]);
  int yt_in       = (int)atoi(argv[10]);
  std::string file_tag  = argv[11];

  // Reference SM Higgs mass for EW precision observables
  double mh_ref = 125.;

  // Create SM and set parameters
  SM sm;
  sm.set_qmass_pole(6, 172.5);
  sm.set_qmass_pole(5, 4.75);
  sm.set_qmass_pole(4, 1.42);
  sm.set_lmass_pole(3, 1.77684);
  sm.set_alpha(1. / 127.934);
  sm.set_alpha0(1. / 137.0359997);
  sm.set_alpha_s(0.119);
  sm.set_MZ(91.15349);
  sm.set_MW(80.36951);
  sm.set_gamma_Z(2.49581);
  sm.set_gamma_W(2.08856);
  sm.set_GF(1.16637E-5);

  // Create 2HDM and set SM parameters
  THDM model;
  model.set_SM(sm);

  bool pset = model.set_param_phys(mh_in,mH_in,mA_in,mC_in,sba_in,l6_in,l7_in,m12_2_in,tanb_in);

  if (!pset) {
    cerr << "The specified parameters are not valid\n";
    return -1;
  }

  // Set Yukawa couplings
  model.set_yukawas_type(yt_in);

  // Prepare to calculate observables
  Constraints constr(model);

  double S, T, U, V, W, X;
  constr.oblique_param(mh_ref, S, T, U, V, W, X);

#if defined HiggsBounds
  HBHS hbhs{};
  const HBHSResult hbhs_result = hbhs.check(model);
//  printf("Hb-obsratio: ",hbhs_result.hb.obsratio[0]);
//  printf("\nConstraints:\n");
//  printf(" Tree-level unitarity: %s\n", constr.check_unitarity());

#endif

  // Prepare to calculate decay widths
  DecayTable table(model);

  //############  h decays #####################################
  double BrhWW = table.get_gamma_hvv(1,3)/table.get_gammatot_h(1);
  double BrhZZ = table.get_gamma_hvv(1,2)/table.get_gammatot_h(1);
  double Brhtautau = table.get_gamma_hll(1,3,3)/table.get_gammatot_h(1);
  double Brhmumu = table.get_gamma_hll(1,2,2)/table.get_gammatot_h(1);
  double Brhss = table.get_gamma_hdd(1,2,2)/table.get_gammatot_h(1);
  double Brhbb = table.get_gamma_hdd(1,3,3)/table.get_gammatot_h(1);
  double Brhuu = table.get_gamma_huu(1,1,1)/table.get_gammatot_h(1);
  double Brhcc = table.get_gamma_huu(1,2,2)/table.get_gammatot_h(1);
  double Brhgg = table.get_gamma_hgg(1)/table.get_gammatot_h(1);
  double Brhgaga = table.get_gamma_hgaga(1)/table.get_gammatot_h(1);
  double BrhgaZ = table.get_gamma_hZga(1)/table.get_gammatot_h(1);
  double BrhAA = table.get_gamma_hhh(1,3,3)/table.get_gammatot_h(1);
  double BrhZA = table.get_gamma_hvh(1,2,3)/table.get_gammatot_h(1);

  // H decays #####################################################
  double BrHWW = table.get_gamma_hvv(2,3)/table.get_gammatot_h(2);
  double BrHZZ = table.get_gamma_hvv(2,2)/table.get_gammatot_h(2);
  double BrHtautau = table.get_gamma_hll(2,3,3)/table.get_gammatot_h(2);
  double BrHbb = table.get_gamma_hdd(2,3,3)/table.get_gammatot_h(2);
  double BrHHpHm = table.get_gamma_hhh(2,4,4)/table.get_gammatot_h(2);
  double BrHhh = table.get_gamma_hhh(2,1,1)/table.get_gammatot_h(2);
  double BrHAA = table.get_gamma_hhh(2,3,3)/table.get_gammatot_h(2);
  double BrHZA = table.get_gamma_hvh(2,2,3)/table.get_gammatot_h(2);
  double BrHWhp =  table.get_gamma_hvh(2,4,4)/table.get_gammatot_h(2);
  double BrHgg = table.get_gamma_hgg(2)/table.get_gammatot_h(2);
  
  // A decays #####################################################
  double BrAtautau = table.get_gamma_hll(3,3,3)/table.get_gammatot_h(3);
  double BrAmumu = table.get_gamma_hll(3,2,2)/table.get_gammatot_h(3);
  double BrAbb = table.get_gamma_hdd(3,3,3)/table.get_gammatot_h(3);
  double BrAZh = table.get_gamma_hvh(3,2,1)/table.get_gammatot_h(3);
  double BrAZH = table.get_gamma_hvh(3,2,2)/table.get_gammatot_h(3);
  double BrAgaga = table.get_gamma_hgaga(3)/table.get_gammatot_h(3);
 
 // hp decays #####################################################  
   double BrHptb = table.get_gamma_hdu(4,3,3)/table.get_gammatot_h(4);
   double BrHpcs = table.get_gamma_hdu(4,2,2)/table.get_gammatot_h(4);
   double BrHpcb = table.get_gamma_hdu(4,3,2)/table.get_gammatot_h(4); 
   double BrHptanu = table.get_gamma_hln(4,3,3)/table.get_gammatot_h(4); 
   double BrHpmunu = table.get_gamma_hln(4,2,2)/table.get_gammatot_h(4);   
   double BrHpenu = table.get_gamma_hln(4,1,1)/table.get_gammatot_h(4);
   double BrHpWh = table.get_gamma_hvh(4,3,1)/table.get_gammatot_h(4);
   double BrHpWH = table.get_gamma_hvh(4,3,2)/table.get_gammatot_h(4);
   double BrHpWA = table.get_gamma_hvh(4,3,3)/table.get_gammatot_h(4);

 // top decays #####################################################
   double BrtHpb = table.get_gamma_uhd(3,4,3)/table.get_gammatot_top();
   double Xs_pp_tHpb = 966.01*BrtHpb;


 // #####################################################  
  std::string filename_param_chisq = file_tag;

  std::ofstream file_param_chisq;
  file_param_chisq.open(filename_param_chisq.c_str(), std::ios_base::app);

  std::stringstream line;
   line << std::setprecision(6) << std::fixed << mh_in << " ";    // 1
   line << std::setprecision(6) << std::fixed << mH_in << " ";    //2
   line << std::setprecision(6) << std::fixed << mA_in << " ";    //3
   line << std::setprecision(6) << std::fixed << mC_in << " ";    //4
   line << std::setprecision(6) << std::fixed << sba_in << " ";   //5
   line << std::setprecision(6) << std::fixed << tanb_in << " ";  //6
   line << std::setprecision(6) << std::fixed << m12_2_in << " "; //7
#if defined HiggsBounds
   line << std::setprecision(6) << std::fixed << hbhs_result.hb.obsratio[0] << " ";//8
   line << std::setprecision(6) << std::fixed <<  hbhs_result.hb.channel[0] << " ";// 9
   line << std::setprecision(6) << std::fixed << hbhs_result.hs.chisq << " ";   // 10 
   line << std::setprecision(6) << std::fixed << hbhs_result.hs.pvalue << " "; //11
#endif
   line << std::setprecision(1) << std::fixed << constr.check_stability() << " "; //12
   line << std::setprecision(1) << std::fixed << constr.check_unitarity() << " ";  //13
   line << std::setprecision(1) << std::fixed << constr.check_perturbativity() << " ";  //14

   //############  h decays #####################################
   line << std::setprecision(6) << std::fixed << BrhWW << " "; //15
   line << std::setprecision(6) << std::fixed << BrhZZ << " ";  //16
   line << std::setprecision(6) << std::fixed << Brhtautau << " "; //17
   line << std::setprecision(6) << std::fixed << Brhmumu << " ";  //18
   line << std::setprecision(6) << std::fixed << Brhss << " "; //19
   line << std::setprecision(6) << std::fixed << Brhbb << " ";  //20
   line << std::setprecision(6) << std::fixed << Brhuu << " "; //21
   line << std::setprecision(6) << std::fixed << Brhcc << " "; //22
   line << std::setprecision(6) << std::fixed << Brhgg << " ";//23
   line << std::setprecision(6) << std::fixed << Brhgaga << " ";//24
   line << std::setprecision(6) << std::fixed << BrhgaZ << " ";//25
   line << std::setprecision(6) << std::fixed << BrhAA<< " "; //26
   line << std::setprecision(6) << std::fixed << table.get_gammatot_h(1) << " "; //27

   //############  H decays #####################################
   line << std::setprecision(6) << std::fixed << BrHWW << " ";   //28
   line << std::setprecision(6) << std::fixed << BrHZZ << " ";   //29
   line << std::setprecision(6) << std::fixed << BrHtautau << " "; //30
   line << std::setprecision(6) << std::fixed << BrHbb << " ";  //31
   line << std::setprecision(6) << std::fixed << BrHHpHm << " ";  //32
   line << std::setprecision(6) << std::fixed << BrHhh << " ";   //33
   line << std::setprecision(6) << std::fixed << BrHAA << " ";   //34
   line << std::setprecision(6) << std::fixed << BrHZA << " ";   //35
   line << std::setprecision(6) << std::fixed << BrHWhp << " ";   //36
   line << std::setprecision(6) << std::fixed << table.get_gammatot_h(2) << " ";  //37

  //##########  A decays  ###########################################
  line << std::setprecision(6) << std::fixed << BrAtautau << " ";  //38
  line << std::setprecision(6) << std::fixed << BrAmumu << " ";   //39   
  line << std::setprecision(6) << std::fixed << BrAbb << " ";  //40     
  line << std::setprecision(6) << std::fixed << BrAZh << " "; //41
  line << std::setprecision(6) << std::fixed << BrAZH << " "; //42
  line << std::setprecision(6) << std::fixed << table.get_gammatot_h(3) << " ";  //43  

  //##########  hp decays  ###########################################
  line << std::setprecision(6) << std::fixed << BrHpmunu << " "; // 44
  line << std::setprecision(6) << std::fixed << BrHpenu << " "; // 45
  line << std::setprecision(6) << std::fixed << BrHpcs << " "; // 46
  line << std::setprecision(6) << std::fixed << BrHpcb << " "; // 47
  line << std::setprecision(6) << std::fixed << BrHpWh << " "; // 48
  line << std::setprecision(6) << std::fixed << BrHpWH << " "; // 49
  line << std::setprecision(6) << std::fixed << BrHpWA << " "; // 50
  line << std::setprecision(6) << std::fixed << BrHptanu << " "; // 51
  line << std::setprecision(6) << std::fixed << BrHptb << " "; // 52

 //##########  t-hp decays  ###########################################
  line << std::setprecision(6) << std::fixed << BrtHpb << " "; // 53
  line << std::setprecision(6) << std::fixed << Xs_pp_tHpb << " "; // 54
  line << std::setprecision(6) << std::fixed << BrHgg << " ";  //55
  line << std::setprecision(6) << std::fixed << BrAgaga << " ";  //56
  line << std::setprecision(6) << std::fixed << S << " "; //57
  line << std::setprecision(6) << std::fixed << T << " ";  //58
  line << std::setprecision(6) << std::fixed << U << "\n";  //59
  file_param_chisq << line.rdbuf();
 
}
